Control Theory Based Airfoil 
Design Using the Euler Equations 


Antony Jameson and James Reuther 


The Research Institute of Advanced Computer Science is operated by Universities Space Research 
Association, The American City Building, Suite 212, Columbia, MD 21044, (410) 730-2656 


Work reported herein was sponsored by NASA under contract NAS 2-13721 between NASA and the 
Universities Space Research Association (USRA). 



Control Theory Based Airfoil Design 
using the Euler Equations 


A. Jameson* 

Department of Mechanical and Aerospace Engineering 
Princeton University 
Princeton, New Jersey 08544, U.S.A. 

and 

J. Reuther * 

Research Institute for Advanced Computer Science 
Mail Stop T20G-5 
NASA Ames Research Center 
Moffett Field, Califonia 94035, U.S.A. 


Abstract 

This paper describes the implementation of op- 
timization techniques based on control theory for 
airfoil design. In our previous work [9, 10, 15] 
it was shown that control theory could be em- 
ployed to devise effective optimization procedures 
for two-dimensional profiles by using the potential 
flow equation with either a conformal mapping or a 
general coordinate system. The goal of our present 
work is to extend the development to treat the Eu- 
ler equations in two-dimensions by procedures that 
can readily be generalized to treat complex shapes 
in three-dimensions. Therefore, we have developed 
methods which can address airfoil design through 
either an analytic mapping or an arbitrary grid per- 
turbation method applied to a finite volume dis- 
cretization of the Euler equations. Here the control 
law serves to provide computationally inexpensive 
gradient information to a standard numerical opti- 
mization method. Results are presented for both the 
inverse problem and drag minimization problem. 

Nomenclature 

Ai,A 2 Cartesian flux Jacobians 
b design variable 

b bandwidth of full Jacobian matrix 
B bounding of flowfield domain at far field 
c, s transformation metric 
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c chord length 

C bounding surface of flowfield domain on airfoil 
C U C 2 transformed flux Jacobians 
Cd coefficient of drag 
Ci coefficient of lift 
C p coefficient of pressure 
Cp coefficient of pressure for sonic flow 
D flowfield domain 
E total energy 

/,<; components of Cartesian fluxes 
f,g /, <7 without pressure terms 
F , G components of transformed fluxes 
F y G F,G without pressure terms 

T control law, physical location of boundary 
Q gradient of design space 
Q projected Q into admissible space 
h modulus of conformal mapping transformation 
H total enthalpy 
I cost function 

J Jacobian of generalized transformation 
K grid transformation matrix 
K f K at the profile surface 
K \ , K 2 grid transformation matrices 

ra number of flowfield evaluations per line search 
M local Mach number 
Moo Mach number at infinity 
n number of design variables 
n number of grid points 
ni ,«2 components of a unit vector normal 
p pressure 
pd desired pressure 
q speed 

q oo speed at infinity 
r radial component in mapped plane 
R generic governing equation for flowfield 
R normalized arc length in t; 
s arc length along airfoil surface 
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5 surface displacement in mapped plane 
t time 

<i,<2 exponents on basis functions 
u , v Cartesian velocity components 
t/, V contravariant velocity components 
w state vector of flowfield unknowns 
x, y Cartesian coordinates 
£, T} body fitted coordinates 
z physical plane 

Z modified boundary condition for V* 
a angle of attack 
ft smoothing parameter 

6 variational 

7 ratio of specific heats 

A distance along search direction 
Ai , A 2 smoothing parameters 

A distance along search direction 
Qi,f2 2 cost function weighting factors 
vector co-state variable 
p density 

poo freestream density 
<x conformally mapped plane 
S angle around circle 


Formulation of the design 
problem as a control problem 

The final efficiency and success of a design de- 
pend on complex multi-disciplinary trade-offs be- 
tween factors such as aerodynamic efficiency, struc- 
tural weight, stability and control, and the volume 
required to contain fuel and payload. Typically a 
design is finalized after numerous iterations, cycling 
between the disciplines. While full multidisciplinary 
optimization is the ultimate goal, a necessary inter- 
mediate step is the development of efficient methods 
for aerodynamic shape optimization. Early investi- 
gations of aerodynamic optimization [16] relied on 
brute force methods, in which the influence of each 
design parameter was estimated by making a small 
change in that parameter and recalculating the flow. 
With a large number of parameters this becomes ex- 
tremely expensive. 

Alternatively, it has been recognized that the de- 
signer has an idea of the kind of pressure distribution 
that will lead to the desired performance. Thus, it is 
useful to consider the inverse problem of calculating 
the shape that will lead to a given pressure distribu- 
tion. This method has the advantage that only one 
flow solution is required to obtain the desired design. 
Unfortunately, unless the pressure distribution satis- 
fies certain constraints a physically realizable shape 
may not necessarily exist. Thus the problem must 
be very carefully formulated. 

The difficulty that the objective may be unattain- 
able can be circumvented by regarding the design 
problem as a control problem in which the control 
is the shape of the boundary. A variety of alterna- 
tive formulations of the design problem can then be 


treated systematically within the framework of the 
mathematical theory for control of systems governed 
by partial differential equations [12]. This approach 
to optimal aerodynamic design was introduced by 
Jameson [9, 10], who examined the design problem 
for compressible flow with shock waves, and devised 
adjoint equations to determine the gradient for both 
potential flow and also flows governed by the Eu- 
ler equations. More recently Ta’asan, Kuruvila, and 
Salas, implemented a one shot approach in which the 
constraint represented by the flow equations need 
only be satisfied by the final converged solution [20]. 
Pironneau has also studied the use of control theory 
for optimum shape design of systems governed by el- 
liptic equations [14], and adjoint methods have also 
been used by Baysal and Eleshaky [2]. 

For flow about an airfoil or wing, the aerodynamic 
properties which define the cost function are func- 
tions of the flow-field variables (u>) and the physical 
location of the boundary, which may be represented 
by the function T , say. Then 


I = l(w,T)> 


and a change in T results in a change 


61 - 


d JL 

dw 


6w + ^6f, 


( 1 ) 


in the cost function. As pointed out by Baysal and 
Eleshaky [2] each term in (1), except for 6w , can 
be easily obtained. J- and can be obtained di- 
rectly without a flowfield evaluation since they are 
partial derivatives. 6T can be determined by either 
working out the exact analytical values from a map- 
ping, or by successive grid generation for each design 
variable, so long as this cost is significantly less then 
the cost of the flow solution. Brute force methods 
evaluate the gradient by making a small change in 
each design variable separately, and then recalculate 
both the grid and flow-field variables. This requires 
a number of additional flow calculations equal to the 
number of design variables. Using control theory, 
the governing equations of the flowfield are intro- 
duced as a constraint in such a way that the final 
expression for the gradient does not require multi- 
ple flow solutions. In order to achieve this 6w must 
be eliminated from (1). The governing equation R 
expresses the dependence of w and T within the 
flowfield domain D , 


R (w, T) — 0, 

Thus 6w is determined from the equation 


6R = 


dR 


dw 


dR 


dT 


6F = 0. 


(2) 


6w + 

Next, introducing a Lagrange Multiplier we have 
61 = 


dF c dF 

—6XV + —6T 


2 



* T ( 

'dR 

dw 

Sw + 

dF_ 

dw 

-ri> T 

'dR " 

dw 


5ft 


(dT_ _ 

\ dT 


r 


d?\ J 


dT 
8w 
8T 


8T 


Choosing to satisfy the adjoint equation 

|T 


5ft 


dw 


w dw 


the first term is eliminated, and we find that 


SI = G J ST 


(3) 

(4) 


where 


g t 



' dR 
dT. ' 


The advantage is that (4) is independent of 6w , with 
the result that the gradient of I with respect to 
an arbitrary number of design variables can be de- 
termined without the need for additional flow-field 
evaluations. The main cost is in solving the ad- 
joint equation (3). In general, the adjoint problem is 
about as complex as a flow solution. If the number 
of design variables is large, the cost differential be- 
tween one adjoint solution and the large number of 
flowfield evaluations required to determine the gra- 
dient by brute force becomes compelling. Instead of 
introducing a Lagrange multiplier, rj), one can solve 
(2) for Sw as 


Sw — — 


'dR' 

-l 

"5ft" 

dw _ 


dT 


(5) 


and insert the result in (1). This is the implicit 
gradient approach, which is essentially equivalent to 
the control theory approach, as has been pointed out 
by Shubin and Frank [18, 19]. In any event, there is 
an advantage in determining the gradient Q by the 
solution of the adjoint equation. 

Once equation (4) is established, an improvement 
can be made by with a shape change 


ST = -A Q 


where A is positive, and small enough that the first 
variation is an accurate estimate of SI . Then 

si = -xg T g < o 

After making such a modification, the gradient can 
be recalculated and the process repeated to follow a 
path of steepest descent until a minimum is reached. 
In order to avoid violating constraints, such as a 
minimum acceptable wing thickness, the gradient 
may be projected into the allowable subspace within 


which the constraints are satisfied. In this way pro- 
cedures can be devised which must necessarily con- 
verge at least to a local minimum. The efficiency 
may be improved by performing line searches to find 
the minimum in a search direction defined by the 
negative gradient, and also by the use of more so- 
phisticated descent methods such as conjugate gra- 
dient or quasi-Newton algorithms. There is the pos- 
sibility of more than one local minimum, but in any 
case the method will lead to an improvement over 
the original design. Furthermore, unlike the tra- 
ditional inverse algorithms, any measure of perfor- 
mance can be used as the cost function. 

Steps (1-4) may be applied to the discrete equa- 
tions which approximate the governing differential 
equations. In the present method they are applied 
instead directly to the differential equations, and the 
resulting adjoint differential equation is then dis- 
cretized in a manner similar to the discretization 
of the flow equations. The former approach is now 
gaining favor in the work of Newman and Taylor el 
al [13, 11] and also Baysal, Eleshaky and Burgreen 
[1, 2, 3, 4]. The two methods can be very similar, 
depending upon the discretization of (3) or (5). The 
current method has an advantage in that the dis- 
cretization and iteration scheme used to solve the 
flowfield system can also be applied directly to the 
adjoint system. Therefore the robust iteration algo- 
rithms and convergence acceleration techniques that 
have been matured for CFD algorithms can be di- 
rectly ported for the solution of the adjoint system. 
Methods applying (1 - 4) to the discrete system of- 
ten resort to matrix inversion methods to solve (3) or 
(5). While direct inversion techniques can be robust 
and reliable they suffer when the number of points 
becomes large because the operation count grows as 
0(n/r), where n is the number of unknowns and 6 
is the bandwidth which can approach 0(n). When 
the number of mesh points becomes large, especially 
in the case of three-dimensional problems, the 0(h) 
operational counts of explicit iteration schemes used 
in CFD can significantly reduce the required time to 
solve the adjoint system. 

In our previous application of control theory to op- 
timal aerodynamic shape design a successful numer- 
ical implementation was demonstrated using the po- 
tential flow equation, with either a conformal map- 
ping or a finite volume discretization. In this work, 
the Euler equations are employed instead of the po- 
tential flow equation to provide a more accurate rep- 
resentation of flows containing strong shock waves 
and entropy variations. Again both an analytic 
mapping and general finite volume discretization are 
explored as possible vehicles for the implementation 
of design via control theory. This is intended to 
be a precursor for the three-dimensional problem, 
where the geometry may be very complex and both 
shocks and vortex roll-up must be considered. The 
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treatment of viscous flows with the Navier Stokes 
equations is a natural future step. 


Design of airfoils using 
the Euler equations 


This section develops the general application of con- 
trol theory for aerodynamic shape optimization us- 
ing the Euler equations. Consider the case for com- 
pressible flow over an airfoil. In the absence of sep- 
aration and other strong viscous effects, the flow is 
well approximated by the Euler equations. In con- 
trast to our previous implementations which relied 
on the isentropic potential equation [15], here strong 
inviscid shocks are modeled correctly with entropy 
production. Consider the flow in a domain D. The 
profile defines the inner boundary C, while the outer 
boundary B is assumed to be distant from the pro- 
file. Let p, p, tt, u, E and H denote the pressure, 
density, Cartesian velocity components, total energy 
and total enthalpy. For a perfect gas 


P= (t- !)p (u 2 + t> 2 ) 

} (6) 

and 


pH = pE + p, 

(7) 

where 7 is the ratio of the specific heats, 
equations may then be written as 

The Euler 

dw df dg . 

(8) 

~dr + di + di = 0 mD ' 


Introduce contravariant velocity components 
= K~ x 


U 

V 


dy dx 

dr} dr] 

dy dx 

di di J 


The Euler equations can be written in divergence 
form as 


with 


dW 8F dG „ 

nr + nr + m =0 mD ' 


W = J 


F = J 


G = J 


P 

pu 
pv 

pE 

pU 

p Uu + Up 

pUv + J£p 

pUH 

pV 

pVu+gp 
pVv + gp 

pVH 


( 10 ) 


( 11 ) 


where x and y are Cartesian coordinates, t is the 
time coordinate and 


w 


f 


P 

pu 

pv ’ 
pE 

pv 
pvu 

pv 2 + p 
pvH 

(9) 


pu 


pu 2 + p 

, 9 = 

puv 


puH 



The coordinate transformations may be defined by 
the transformation matrix 

dx dx 
di Ef) 
dxt 

di dfj _ 



with £ and r; representing the computational coor- 
dinate system and the Jacobian 


r , / r ^ \ 9x dy 


dx dy 
dr] di 


Assume now that the computational coordinate sys- 
tem conforms to the airfoil section in such a way 
that the surface C is represented by p = 0. Then 
the flow is determined as the steady state solution 
of the equation (11) subject to the flow tangency 
condition 

V = 0 on C. (12) 

At the far field boundary B , conditions are speci- 
fied for incoming waves, while outgoing waves are 
determined by the solution. 

As a first example of the use of control theory, 
consider the case of the inverse problem where the 
cost function may be defined as 

; = i/ c ( P -P. )! * = i/ c ( P -p 1 ) J 

where p<* is the desired pressure. The design problem 
is now treated as a control problem where the control 
function is the airfoil shape, which is to be chosen 
to minimize I subject to the constraints defined by 
the flow equations (10-12). A variation in the shape 
will cause a variation Sp in the pressure in addition 
to a variation in the geometry and consequently the 
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variation in the cost function becomes 


On the profile n\ — 0 and n 2 = — 1. It follows from 
equation (12) that 


81 = 


+ 


J^{p-Pd)6p 

< 13 > 


Since p depends on w through the equation of 
state (8-9), the variation Sp can be determined from 
the variation Sw. Define the Jacobian matrices 


A\ 


d£ 

dw ' 


A 2 


dg_ 
dw 1 


Ci = '£, J K^Aj. (14) 
j 


Then the equation for Sw in the steady state be- 
comes 

±(6F) + -?-(6G) = 0, (15) 

where 

6F = C,6w + 6^J^jf + 6^J^jg 

SG = c^ + *(^)/ + «(j|) 5 . 

Now, multiplying by a vector costate variable ip and 
integrating over the domain we have 




d^drj = 0. 


0 


0 


+ p 

‘(jS) 

f | ] 6 p 


0 


0 


(16) 


Suppose now that ip 
of the adjoint equation 


is the steady state solution 


dip 

dt 


r ^T dip n t 9lp 

Gl 


= o 


in D. 


(17) 


At the outer boundary incoming characteristics for 
ip correspond to outgoing characteristics for Sw. 
Consequently, one can choose boundary conditions 
for ip such that 


riiip T CiSw — 0 . 


Then if the coordinate transformation is such that 
S ( J 1 ) is negligible in the far field, the only re- 
maining boundary term is 

L* tsg %« 

Thus by letting ip satisfy the boundary condition, 


If ip is differentiable this may be integrated by parts 
to give 

L{% iF+ % 6G ) d ^ = 

/ (n x t P T 6F + h 2 ip T 6G)d(, 

Jb 

+ J (n i ij> T 6F + n 2 rJ> T 6G)dZ, 

where n : are components of a unit vector normal to 
the boundary. No boundary integrals appear in the 
r } direction because the mesh is assumed to be of O- 
type, with the result that the solution is periodic in 
the £ coordinate thereby cancelling the 7] boundary 
integrals. Thus the variation in the cost function 
may now be written 


SI — 


+ 

+ 


J c (p-Pd)bp d£ 

[ ( ^ >T 

Jd v ae 


T dii> T 

6F+-?—6G)dtd v 


[ (ii^ t 6F + n 2 ip T 8G)d{, 

Jb 

I (n\ip T 6F + n 2 t/> T 6G ) d£. 

Jc 


on C„ (18) 

we find finally that 


SI = 


+ 


+ 


-f 


-f 


\ J {p-Pd) 2 S d{ 

fMmw m 

l S c (r ■ rdf s (I) * 

f ( x ( ® y \ t x ( dx \ ^ 

Jd \ V# 1 // \dTj) 9 ) 

♦£K 8 ) /+ ‘ (*)')«* 
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If the flow is subsonic, this procedure should con- 
verge toward the desired pressure distribution since 
the solution will remain smooth, and no unbounded 



derivatives will appear. If, however, the flow is tran- 
sonic, one must allow for the appearance of shock 
waves in the trial solutions, even if pa is smooth. In 
such instances p — p<* is not differentiable. This dif- 
ficulty can be circumvented by a more sophisticated 
choice of the cost function. Consider the choice 


the derivative of the mapping function be 

£ = he‘». 
da 

Now using polar coordinates r, and 6 in the a plane, 
the first transformation matrix is 



where A] and A 2 are parameters, and the periodic 
function Z(£) satisfies the equation 


X\Z 



= P-Pd> 


( 20 ) 


Then, 

" = l c { x ' zsz+ ^H 6z )%^ 

= S C Z {’" 6Z -^ 6Z )%^ 

Thus, Z replaces p — pd with a corresponding mod- 
ification to the boundary condition for the adjoint 
equations. 

The following two sections develop alternative im- 
plementations of the general method discussed thus 
far. The differences between the two lie in the way in 
which the control modifies the mesh on which both 
the flow and the adjoint equations are solved. The 
first method relies on an analytic mapping of the 
mesh and directly uses the value of the mapping at 
the surface as the control. This implies that each 
surface mesh point becomes a design variable. The 
second method allows for an arbitrary numerically 
generated initial mesh. This is used as a template 
for all subsequent meshes, which are created by a 
grid perturbation technique. A set of design vari- 
ables spanning a space of allowable airfoils is then 
chosen to assure smooth variations in the geometry. 
Thus these design variables become the control, and 
no dependence upon transformations is necessary. 


Design Using an Analytic Mapping 


A convenient way to treat an airfoil is to use a con- 
formal mapping of the profile in the z plane to a 
near circle in the a plane, followed by shearing of 
the radial coordinate to make the system bound- 
ary conforming. Polar coordinates are introduced 
in the mapped plane a . When mapped back to the 
physical plane this gives a smooth, nearly orthog- 
onal grid. We can now specialize our generalized 
design procedure to treat this grid system. Define 
the first conformal mapping from z to a by letting 


*0 

X r 

= h 

rs 

c 

. v* 

Vr 


— cs 

s 


K x = 


and we can define contravariant velocities 


U 

V 


s — c 

c s 


][: 


where 


s — sin (/? — 0 ) , c — cos (/? — 0) . 


The Euler equations can now be represented in the 
a plane as 


d(rh*W) . d(hF) , d(rhG) n . n 
+ — ™ 1- — = 0 in D, 


dt 


89 


dr 


( 21 ) 


where 


W = 


P 

pu 

pv 

pE 


pU 


P V 

pU u 4 sp 

. G = 

pV u 4 cp 

pU v — cp 


pV v 4 sp 

pUH 


pVH 


F = 


( 22 ) 

Now let the final computational coordinates be de- 
fined by a radial shearing transformation 

* = r = 9 + S(0 

and the transformation matrix 


Ko = 


r de 

ae 


1 

0 

FI 




dr 

dr 


dS 

1 


d V _ 




Ji = T 


Now we can identify the complete transformation 
matrix as 


K = K\ K 2 = h 
while the fluxes are 


rs + S{C c 
[ -rc -f SfS s J ’ 


hF = h(sf-cg) 
~ Vt}/ ~ %r\9 
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and 


h[(t} + S)c- S ( s]f 
+h [(*? + <5) s + S( c] g 
- Vif. 


Then 


61 


h[(ij + S)G-S ( F] = 


Thus the Euler equations assume the form 


= J (Q6S-P6S ( )dt 


- L 

where the gradient is 


G6Sd£, 


1L(( V + S)h>w) 

+ I<W) 

+ -jt(h(ti + S)G-hS ( F) = 0, 

while the surface tangency condition on the velocity 
becomes 

X{V — y^u = h [(;; + s) V — s ( U] = 0. 

Now we take S(£) as the control. It is also conve- 
nient to represent the inverse problem by the cost 
function 

/ = ^ J (p - p d f d0 = /= ^ J^ip- Pdf d£. 

This eliminates terms in 6 from the gradient. 
The variations in the fluxes become 


6{hF) = C\6w 

6[h{ri + S)G- hS^F] = C 2 Sw + h6SG - hSS^F 

where C 1 and C 2 are the Jacobian matrices defined 
in equation (14). Choosing V' to satisfy the adjoint 
equation (17) with the boundary condition 


X£tl > 2 - 3 = h [(r) + S) s + S^c] tl>n = P-Pd 

the variation in the cost reduces to 

61 = f (P~ Prf ) 6pd Z 

+ jj T {k F+ k G ) d(d " 


f iP T (SShG - 6S{hF) d£ 

Jc 

+ L 


r\ 

rP T ( 6S hG + 6S ( hF) didr, 

oi 


where F and G are the fluxes defined in equation 
(22), and F and G are F and G with the pressure 
terms deleted. Define 


P — ij? hF + J tp 7 ^(hF)drf 


d_ 

dr} 


Q = i> T hG + J \p T -j^(hG)di,. 


(23) 


G = Q+ ~di' (24) 

Here, the gradient is defined on the surface of the 
computational mesh. In effect, the gradient with 
respect to a number of design variables equal to the 
number of surface points can be calculated in one 
shot. 


Design Using an Arbitrary Mesh 


In order to treat a more general mesh we revert 
to equations (13-19). The difficulty in using these 
general equations is that the variation of the met- 
ric terms in the equations needs to be obtained in 
order to construct 61 in equation (19). One way 
to accomplish this is by using an existing grid gen- 
eration algorithm and applying finite differences to 
calculate the necessary information. While this ap- 
proach would still obviate the use of multiple flow 
solutions to determine the gradient, it would still re- 
quire the mesh generator to be used repeatedly. The 
number of mesh solutions required would be pro- 
portional to the number of design variables. This 
may be acceptable, since the flow solution process 
is typically much more computationally expensive 
then grid generation. Such a method should then 
ensure a significant savings over using finite differ- 
ences for both the grid generation and flow solu- 
tion processes, and it was successfully applied by 
the authors to two-dimensional problems in poten- 
tial flow [15]. For three-dimensional design, where 
both the number of design variables and the com- 
putational cost of grid generation can be high, the 
method is excessively expensive. Further, for com- 
plicated three-dimensional configurations, for which 
fully automatic grid generation is still not practical, 
it will not be feasible. 

This motivates the need to find a method which 
bypasses these difficulties. In order to remove 
the cost of the successive grid generation from the 
gradient calculation a successive grid perturbation 
method is therefore used. In this approach, which 
was also used by Burgreen and Baysal [4], an initial 
structured curvilinear body fitted grid over the ini- 
tial configuration is created by any grid generation 
process before optimization. Then the geometry as 
well as the grid become inputs to the optimization 
process. New grids, which conform to the surface 
as it is modified, can then be generated by shifting 
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the grid points along each grid index line projecting 
from the surface by an amount which is attenuated 
as the arc length from the surface increases. If the 
outer boundary of the grid domain is held constant 
the modification to the grid has the form 

x neu, _ x ol d + n ( x new _ 

l r w = y° ld +n(y?™-y°“), 

where x and y represent the volume grid points and, 
x s and y s represent the surface grid points. 72 rep- 
resents the arclength along the surface projecting 
mesh line measured on the original grid, from the 
outer domain, and normalized so that 72 = 1 at the 
inner surface. The required variations in the met- 
ric terms can then be obtained in terms of surface 
perturbations since it follows that, 


6x = 7 Z6x s 

Sy = 7 Z6y s . 


and 

[K] = n[K 9 ] (25) 

where K s is the transformation matrix K defined at 
the surface. 

Now it is convenient to rewrite (19) after integrat- 
ing by parts as. 

"= Ua (P ~ Pi)2S (^) d( 



(26) 


where / and g are again the flux components / and 
g with the pressure terms deleted from the momen- 
tum equation. Substituting the expressions defined 
by equation (25) into (26) allows us to integrate 
the second term along the index direction projecting 
from the configuration, r) surface without any depen- 
dence on particular design variables since the metric 
variations are fully determined by the surface per- 
turbations. Thus, the expression for the variation in 
the cost function can be reduced to surface integrals 
only 


61 = 


^( P -Pd) 2 «( 0)* 

f dx\ 
1 “ M "af 


I s (t) 

Jc \ d vJ , 
c ( dy\ 

Jc* T {- S (j%) f + Sl 


Mj 


+ -V 22 dS 


fdx 


) # } 


4, 


( 27 ) 


where 

f ip T -§£ (Rf) di] (Kg) dr) 

f i’ T ^ (72/) dr] f ip T £ (Kg) drj _ 

(28) 

While this type of grid perturbation method does 
not guarantee that grid lines will not eventually 
cross if the perturbations are large, this point is 
irrelevant for gradient calculations since only ana- 
lytic grid derivative information is needed. However, 
since we employ a numerical optimization algorithm 
with line searches along a descent direction, true re- 
grid ding is also necessary. For these line search cal- 
culations the grid perturbation algorithm is used so 
long as negative volumes are not created. If singu- 
larities begin to develop in the grid, the original grid 
generator can be used to create a new grid and the 
process restarted. In practice, the grid perturbation 
method has proven to be robust, and no failures due 
to singularities have occured during optimization of 
typical convex geometries. 

Equation (27) has reduced the expression for 61 
into line integrals along the surface where the only 
remaining unknowns are the grid metrics. These 
surface grid metrics can be easily determined for 
any modification in the surface by direct evalua- 
tion. This suggests choosing a set of design variables 
which smoothly modifies the original shape, say 6;. 
The gradient can then be defined with respect to 
these design variables as 

<?(&<) = ^, (29) 

where 61 is calculated by (27) with each term 6* 
being independently perturbed by a finite step. 

Therefore, to construct Q , an independent basis 
space of perturbation functions fc t , i = l,2,...,n 
(n = number of design variables) is chosen that al- 
lows for the needed freedom of the design space. De- 
sign variables are chosen with the following form, 
suggested by Hicks and Henne [6, 7]: 

b(x) = sin ^ 7 xx io *io(M 
b(x) = x t '(l-x)e~ t > x ) 

where t\ and <2 control the center and width of the 
perturbation and x is the normalized chord length. 
When distributed over the entire chord on both up- 
per and lower surfaces these analytic perturbation 
functions admit a large possible design space. They 
have the advantage of being space based functions, 
as opposed to frequency based functions, and thus 
they allow for local control of the design. They can 
be chosen such that symmetry, thickness, or volume 
can be explicitly constrained. Further, particular 
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choices of these variables will concentrate the design 
effort in regions where refinement is needed, while 
leaving the rest of the airfoil section virtually undis- 
turbed. The disadvantage of these functions is that 
they do not form a complete basis space, nor are 
they orthogonal. Thus, they do not guarantee that 
a solution, for example, of the inverse problem for a 
realizable target pressure distribution will necessar- 
ily be attained. Here, they are employed due to their 
ease of use and ability to produce a wide variation 
of shapes with a limited number of design variables. 

Implementation of the Euler based 
design methods 

Both of the design methods have been successfully 
implemented. The two techniques share many com- 
mon features such as the flow and adjoint solution 
algorithms. The procedures can be summarized as 
follows. 

1. Solve the flow equations (8-12) for />, tz, v, p, 
£, H, U y and V. 

2. Smooth the cost function if necessary by (20). 

3. Solve the adjoint equations (17) for xf) subject 
to the boundary condition (18). 

4. Either calculate P and <5, by (23), from the 
variation in the control S (£), or evaluate Afij 
from equation (28). 

5. Evaluate Q by equation (24) or (29). 

6. Project Q into a feasible direction subject to any 
geometric constraints to obtain Q. 

7. Then either correct the mapping 5(£) or update 
the design variables 6, based on the direction 
from steepest decent 


non-linear aerodynamic design of complete aircraft 
can be brought into the realm of computational fea- 
sibility. The method also has the advantage of be- 
ing quite general in that arbitrary choices for both 
the design variables and optimization technique are 
admitted, as is demonstrated by the two implemen- 
tations. 

The practical implementation of the design meth- 
ods rely heavily upon fast and accurate solvers for 
both the state (iu) and co-state (rft) fields. Further, 
to improve the speed and realizability of the meth- 
ods, a robust choice of the optimization algorithm 
must be made. In this work, Jameson’s FLOS2 
computer program has been used to solve the Eu- 
ler equations. This program uses a multistage time 
stepping scheme with multigrid acceleration to ob- 
tain very rapid steady state solutions, typically in 
25 steps. The adjoint equations are solved by a sim- 
ilar method, in which the flux calculations for the 
Euler equations are replaced by the corresponding 
formulas for the adjoint equation. 

In the implementation with analytic mapping a 
simplified gradient procedure is used as the opti- 
mization process. To preserve the smoothness of the 
profile the gradient is smoothed at each step. Thus 
the change in the shape function S (£) is defined by 
solving 

tS-%zP%z6S = -\G, 
ot oZ 

where ft is a smoothing parameter. Then, to first 
order, the variation in the cost is 


61 = 


< 



di 

d(. 


0. 


6S ( 0 = -A g or Sbi = -A£ 
or as an alternative a quasi-Newton method. 

8. Return to 1. 

In practice these methods resemble those used by 
Hicks, Reuther et al. [16, 8, 17] with control the- 
ory replacing the brute force, finite difference based, 
gradient calculation. In their most recent applica- 
tions, single point design is accomplished on realistic 
wing body configurations with full nonlinear model- 
ing of supersonic flow. Unlike these procedures, the 
current methods’ computational cost do not hinge 
upon the number of design variables, which in these 
cases is either the number of surface mesh points 
used to represent S (f), or the number of perturba- 
tion functions 6*. Thus it can be hoped that with the 
implementation of this method in three-dimensions, 


For the implementation on arbitrary meshes a quasi- 
Newton optimization method is employed. For this 
purpose the QNMDIF program developed by Gill, 
Murray and Wright [5] is used. 

The option to minimize the pressure drag coeffi- 
cient is realized in both methods by redefining the 
cost function as 


I — Cd — 


1 f 

1 2 - / 

i pooqlocJc 


where c is the chord length. To prevent the pro- 
cedure from trying to reduce drag by reducing the 
profile to a non-lifting flat plate a target pressure dis- 
tribution can be retained in the cost function, which 
becomes 

i = \niJ c ( P -p d ) 2 dt + n 2 c d , 

where ft, and Q 2 are weighting parameters. 
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Numerical tests of the design method 
Using an Analytic Mapping 

Three preliminary test cases are presented for the 
design algorithm using an analytic mapping. The 
first two address the problem of attaining a desired 
pressure distribution for lifting airfoils. The most 
convenient method of obtaining such solutions with 
the present design method is to determine the lift 
coefficient associated with the target pressure dis- 
tribution, and match this lift with the initial air- 
foil. Therefore the angle of attack is adjusted after 
every 5 iterations of the flow calculation to drive 
the solution toward the target lift coefficient. The 
first example using this technique, shown in Fig- 
ure 1, drives the Korn airfoil toward the target pres- 
sure distribution for the NACA 64A410 airfoil at 
Moo — 0.75, a = 0°, and C\ — 0.66. This case re- 
quires a shift in the shock location and a significant 
change in the profile shape. The final solution al- 
most exactly recovers the pressure distribution and 
the airfoil shape. In the next example, Figure 2, 
the process is reversed to modify a NACA 64A410 
airfoil to recover the Korn pressure distribution at 
Mach 0.75 as a target. Again the target is almost 
attained, except at the trailing edge, where a per- 
fect cusp is not easy to reproduce. The last test 
case introduces drag as the cost function. Again the 
design process is carried out in the fixed lift mode. 
In Figure 3, a NACA 64A410 is again used as a 
starting airfoil. The design is at M 0 Q = 0.75 and 
C\ = 0.68 where a strong shock causes considerable 
wave drag in the initial airfoil. The objective is to 
reduce the drag with the smallest possible change 
in the airfoil. Therefore the pressure distribution of 
the initial airfoil is retained as the target pressure 
distribution, and the cost function is a blend of the 
drag coefficient and the deviation from the target 
pressure. The final design has a reduction in drag 
from C d = 0.0127 to C d - 0.0016. 


Numerical tests of the design method 
Using a General Mesh 

Several test cases are conducted with the design 
method using the general mesh implementation. 
This technique more closely resembles brute force 
optimization methods. First, perturbations are 
made to individual design variables, and second, the 
calculated gradient vector is fed into a quasi-Newton 
optimization method. This similarity allows for easy 
comparisons of efficiency and accuracy between the 
brute force method and the control theory method. 
A fair test of the efficiency is conducted by calculat- 
ing the gradient using both methods while adjusting 


the level of convergence of the flowfield and adjoint 
solvers to obtain the optimum efficiency for each 
method. For the brute force calculation, the flow 
solver is restarted from a previously converged solu- 
tion and thus requires only a few additional multi- 
grid iterations for each component of the gradient. 
Nevertheless, a dramatic benefit in efficiency is seen 
by using control theory to calculate the gradient. 
Using 50 design variables, 96.5 seconds of Cray C90 
time is required to calculate the brute force gradient, 
compared to only 14.1 seconds required for the cor- 
responding calculation using control theory. With 
the present implementation, the control theory pro- 
duces a less accurate estimate of the gradient than 
the brute force method. Despite this difficulty the 
estimate of the gradient is still more than adequate 
to drive the numerical optimization procedure to- 
ward convergence. 

The first test case for the design method with an 
arbitrary mesh is shown in Figure 4. The NACA 
0012 airfoil is is modified to achieve the pressure dis- 
tribution of a NACA 64A410 airfoil at Mach 0.735, 
a = 0° and C\ — 0.6623. 50 Hicks-Henne design 
variables distributed around the airfoil are used to 
modify the shape with the final airfoil and pressure 
distribution matching the target very closely. The 
second test case also uses 50 design variables and the 
NACA 0012 as a starting condition, as illustrated in 
Figure 5. A target pressure distribution of the RAE 
2822 airfoil at a = 1°, Mach 0.75 and C\ = 0.6982 is 
specified. Again the design method converges to the 
target almost exactly matching even the shock posi- 
tion and strength. A third test case of the method in 
inverse mode is displayed in Figure 6. This time, the 
NACA 0012 airfoil is driven towards the target pres- 
sure distribution of the GAW 72 airfoil operating at 
Mach 0.70, a = 1° and C\ — 0.8158. Under these 
conditions the target airfoil displays a very strong 
shock, yet the design method is able to converge to 
the desired shape without visible discrepancies using 
50 design variables. The fourth test case represents 
a greater challenge. The Korn airfoil at Mach 0.75, 
o = 0° and C\ — 0.6249 is chosen as the target pres- 
sure. The challenge is presented by the fact that the 
Korn airfoil at this condition has shock free super- 
critical flow with no wave drag. Figure 7 shows the 
NACA 0012 airfoil being driven towards the desired 
target, with the design method employing 52 design 
variables. As the design is approached there is a ten- 
dency to produce a double shock pattern instead of 
a smooth recompression. In effect, the design space 
for this problem is more nonlinear than in the pre- 
vious test cases. The final design is very close, but 
does show some discrepences from the desired pres- 
sure distribution. In the fifth test case the RAE 2822 
airfoil is revisited for a target pressure distribution. 
However, a potential flow solver was used to provide 
the target pressure distribution. Thus this pressure 
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distribution is not realizable by the Euler equations 
because the shock wave of the target is modeled as 
an isentropic jump. The final solution with 50 de- 
sign variables, shown in Figure 8, very closely ap- 
proximates the desired pressure distribution, but of 
course does not match it exactly. An examination 
of the final airfoil reveals a striking difference be- 
tween it and the expected airfoil. It can be seen 
that for such strong shock cases, the potential flow 
equation can give quite incorrect results. The sixth 
test case introduces drag as the objective function to 
be minimized. The NACA 64A410 airfoil at Mach 
0.75, a = 1° Ci = 0.700 and Cd = 0.0162 is used as 
a starting condition. Figure 9 illustrates by choos- 
ing 18 design variables that modify the camber, the 
optimization procedure is able to reduce the drag to 
just 7 counts in 3 design cycles. This is accomplished 
while the lift coefficient, Mach number and thickness 
distribution remain unchanged. In the final test case 
drag is again used as the objective function except 
that it is augmented by a lower limit on minimum 
allowable C p . This constraint becomes an additional 
term in the boundary condition of the adjoint equa- 
tion and demonstrates the versatility of the method. 
With this constraint added in the way of a penalty 
function, Figure 10 shows that the NACA 0012 air- 
foil operating at Mach 0.75 and C\ = 0.6 can be 
substantially improved. Twenty design variables are 
chosen which preserve symmetry in the airfoil. With 
a limit of - 1 .2 set for C p , the design procedure is able 
to reduce the drag from 266 counts to 54 counts in 4 
design cycles. The final pressure distribution is ob- 
served to display a flat top at C p = — 1.2 as would 
be expected. 

Conclusions and Recommendations 

We have developed two control theory based airfoil 
design methods for the Euler equations. The differ- 
ence in the two methods lie in the manner in which 
the variations with respect to the position of the 
mesh points are treated. The first relies on an an- 
alytic mapping, while the second permits an arbi- 
trary mesh. The methods represent an extension of 
our previous work on the design of airfoils with the 
potential flow equation. The new methods are both 
efficient and robust, combining the versatility of nu- 
merical optimization methods with the efficiency of 
inverse design methods. A demonstrated seven fold 
improvement in efficiency was realized over the brute 
force optimization method. In three-dimensional 
problems for which the number of design variables 
must be much larger, the improvement in efficiency 
will be even more significant. The direct extendibil- 
ity of this method to three-dimensions is perhaps its 
most appealing aspect, with the final goal being to 
create practical aerodynamic shape design methods 
for complete aircraft configurations. 
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la: Initial Condition 

C { = 0.7019, C d = 0.0015, a = 0.266° 



lb: 50 Design Iterations 

Ci - 0.6610, C d = 0.0136, a = -0.039° 


Figure 1: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
— , x Initial Airfoil: Korn. 

- - + Target C p : NACA 64A410, M = 0.75. 

Inverse Design 



Figure 2: Lifting Design Case, M — 0.75 Fixed Lift Mode 
— , x Initial Airfoil: NACA 64A410. 

, + Target C p : Korn, M = 0.75. 

Inverse Design 
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3a: Initial Condition 3b: 25 Design Iterations 

Ci = 0.6774, C d = 0.0144, a = -0.093° C/ = 0.6855, = 0.0010, a = -0.722° 


Figure 3: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
— , Initial Airfoil: NACA 64A410. 

Drag Reduction 



4a: Initial Condition 

Ci = 0.6623, C d = 0.0281, a = 2.853° 



4b: 22 Design Iterations 

Ci = 0.6623, C d = 0.0079, a = 0.059° 


Figure 4: Lifting Design Case, M = 0.735, Fixed Lift Mode. 
— , x Initial Airfoil: NACA 0012. 

, + Target C p : NACA 64A410, M = 0.735. 

Inverse Design 
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5a: Initial Condition 

Ci = 0.6982, C d = 0.0379, a = 2.593° 



Figure 5: Lifting Design Case, M — 0.75, Fixed Lift Mode. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p : RAE 2822, M = 0.75. 

Inverse Design 



6a: Initial Condition 

Ci = 0.8158, C d = 0.0348, a = 3.835° 
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6b: 32 Design Iterations 
Ci = 0.8158, C d = 0.0159, a = 0.094° 


Figure 6: Lifting Design Case, M = 0.70 r Fixed Lift Mode. 
— , x Initial Airfoil: NACA 0012. 

- - -, + Target C p GAW 72, M = 0.70. 

Inverse Design 
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7a: Initial Condition 

Ci = 0.6249, C d = 0.0294, a = 1.958° 


7b: 21 Design Iterations 

Ci = 0.6249, = 0.0004, a = 0.248° 


Figure 7: Lifting Design Case, Af = 0.75, Fixed Lift Mode. 
— , x Initial Airfoil: NACA 0012. 

, + Target C p : Korn, M — 0.75. 

Inverse Design 
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8b: 21 Design Iterations 

Ci = 0.7970, Cd = 0.0147, a = 1.126° 


Figure 8: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
— , x Initial Airfoil: NACA 0012. 

, + Target C p : RAE 2822 (Potential Flow), M = 0.75. 

Inverse Design 
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9a: Initial Condition 

Ci = 0.7000, C d = 0.0162, a = 0.361° 




9b: 3 Design Iterations 

Ci = 0.7000, C d = 0.0007, a = 0.964° 


Figure 9: Lifting Design Case, M = 0.75, Fixed Lift Mode. 
— , Initial Airfoil: NACA 64A410. 

Drag Reduction 




10a: Initial Condition 

Ci = 0.6000, Cd = 0.0266, a = 2.629° 



10b: 4 Design Iterations 

C) = 0.6000, C d = 0.0054, a = 2.713° 


Figure 10: Lifting Design Case, M — 0.75, Fixed Lift Mode. 
Initial Airfoil: NACA 0012. 

Drag Reduction 
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